# include "scene.h"
# include <stdio.h>
# include <stdlib.h>
# include <math.h>

/* Total number of grid-points. */
#define NXPT 15
/* NDERV-th derivative. */
#define nDERV 2
/* Discretization order of the NDERV-th derivative. */
#define nDISCR 2

extern short cmtrx_nthderv_CD(double *x, unsigned short nxpt, double *xderv,\
 unsigned short ndiscr, struct cmtrx *cm, unsigned short nderv, double **tempmtrx[3]);


int main()
 {
  short nxpt = NXPT;
  short i,j;
  double x[NXPT], u[NXPT], derv, *xderv2, xderv[NXPT-1], temp;
  struct cmtrx cm;

  for(i = 0; i < nxpt; i++)
   {
    x[i] = asin(0.91*cos(2*acos(0)*(nxpt-i)/(double)nxpt))/asin(0.91);
    u[i] = exp(-x[i]);
    if(i > 0) xderv[i-1] = 0.5*(x[i]+x[i-1]);
   }

  xderv2 = 0;

  i = cmtrx_nthderv_CD(x, nxpt, xderv2, nDISCR, &cm, nDERV, (double ***)0);
  printf("ERROR Code = %d\n", i);

//    for(i = 0; i < cm.nrows; i++) { for(j = 0; j < nxpt; j++) printf("%E\t", cm.mtrx[i][j]); printf("\n");}

  if(xderv2 == 0) xderv2 = x;
  for(i = 0; i < cm.nrows; i++)
      {
       for(derv = 0, j = 0; j <= cm.iend[i]-cm.istart[i]; j++) derv += cm.mtrx[i][j]*u[cm.istart[i]+j];
//        for(derv = 0, j = 0; j < nxpt; j++) derv += cm.mtrx[i][j]*u[j];
       temp = 0;
       if((i > 0) && (i < cm.nrows-1)) temp = 2*(((u[i+1]-u[i])/(x[i+1]-x[i]))-((u[i]-u[i-1])/(x[i]-x[i-1])))/(x[i+1]-x[i-1]);
       printf("%0.15E\n%0.15E : du/dx[n=%d]\n%0.15E : Real Value\n\n", temp, derv, i, pow(-1,nDERV)*exp(-xderv2[i]));
      }

  for(i = 0; i < cm.nrows; i++) free(cm.mtrx[i]); free(cm.mtrx);
  free(cm.istart); free(cm.iend);

  return 1;
 }
